library(knitr)
knitr::opts_chunk$set(cache = TRUE, cache.lazy = FALSE, warning=FALSE,
fig.width=6, fig.height=6, autodep=TRUE)
buildFault = function(l=400, setoff=0, rundir="TPV14", run=1) {
# through fault
x0 = 0
y0 = 0
z0 = 0
L = 28e3
W = 15e3
w = l
nl = round(L/l)
nw = round(W/w)
strike = 0
dip = 90
rake = 180
ddot = 0
sectNum = 1
system2("/Users/dinger/RSQSim.github.master/BuildFaultModels/singlePlanarFault",
args=c(x0, y0, z0, L, W, nl, nw, strike, dip, rake, ddot, sectNum),
stdout=file.path(rundir, "fault1.flt"))
# branch fault
strike = 30
x0 = setoff * sin(strike*pi/180)
y0 = 16e3 + setoff * cos(strike*pi/180)
z0 = 0
L = 12e3
nl = round(L/l)
nw = round(W/w)
strike = 30
sectNum = 2
system2("/Users/dinger/RSQSim.github.master/BuildFaultModels/singlePlanarFault",
args=c(x0, y0, z0, L, W, nl, nw, strike, dip, rake, ddot, sectNum),
stdout=file.path(rundir, "fault2.flt"))
system2("cat", args=c(file.path(rundir, "fault1.flt"),
file.path(rundir, "fault2.flt")),
stdout=file.path(rundir, paste0("faults.flt")))
unlink(c(file.path(rundir, "fault1.flt"),
file.path(rundir, "fault2.flt")))
flt = readFault(file.path(rundir, "faults.flt"))
return(flt)
}
Note that the 3 km x 3 km nucleation region centered on \(y\) = 8 km, \(z\) = -7.5 km may or may not fit exactly in a mesh with a given side length.
initStresses = function(flt, sigma=rep(120, length(unique(flt$iseg))),
tau=rep(70, length(unique(flt$iseg))),
ynuc=8000, znuc=-7500, wnuc=1500, segnuc=1, taunuc=81.6) {
sigma.out = sigma[flt$iseg]
write(sigma.out, file=file.path(rundir, "initSigma.txt"), ncol=1)
tau.out = tau[flt$iseg]
nuc = which(flt$iseg == segnuc & abs(flt$z - znuc) <= wnuc & abs(flt$y - ynuc) <= wnuc)
tau.out[nuc] = taunuc
write(tau.out, file=file.path(rundir, "initTau.txt"), ncol=1)
}
These are simple and the same for all the runs we’ll be doing here.
TPV14 (and TPV15) use \(\rho = 2670 \mathrm{kg}/\mathrm{m}^3\), \(V_s\) = 3464 m/s, and \(V_p\) = 6000 m/s. RSQSim needs the Lame parameters \(\lambda\) and \(\mu\). In terms of the density and wave speeds these are
\[\begin{align} \mu &= \rho V_s^2 \\ \lambda &= \rho \left( V_p^2 - 2V_s^2 \right) \end{align}\]rho = 2670
Vs = 3464
Vp = 6000
lameMu = rho * Vs*Vs / 1e6 # in MPa
lameLambda = rho * (Vp*Vp - 2*Vs*Vs) / 1e6
We will use our usual fixed values of \(V^\mathrm{EQ}\) (though will definitely want to try variable \(V^\mathrm{EQ}\) eventually) of 1 m/s (update: given the specified \(V_s\), \(\Delta\tau\), and shear modulus, the shear impedance relation implies a slip speed of 1.57 m/s, but I have adjusted it to 1.2 m/s to better match the rupture propagation speed seen in the dynamic models), \(V^\ast\) of \(10^{-6}\) m/s, \(D_c\) of \(10^{-5}\) m, \(a\) of 0.01, and \(b\) of 0.015. Then we will choose \(\mu_0\), and \(\theta_0\) from the two equations in http://rsqsim.ucr.edu/wiki/index.php/Stress_drops_in_artificially_nucleated_events#Further_attempt_at_more_accurate_estimate_of:
\[\begin{align} \mu_0 &= \frac{\tau_f}{\sigma} + \left(b-a\right)\ln\left(\frac{V^\mathrm{EQ}}{V^\ast}\right) \\ \theta_0 &= \left(\frac{D_c}{V^\ast}\right)^\frac{b-a}{b-a+a^\prime} \left(\frac{l b \sigma}{G V^\mathrm{EQ}}\right)^\frac{a^\prime}{b-a+a^\prime} \exp\left(\frac{\tau_p/\sigma - \mu_0}{b -a +a^\prime}\right) \end{align}\]And then force the nucleation region to fail immediately at its initial shear stress level by lowering \(\theta\) there by the appropriate amount. RSQSim can do that automatically, but I’ll just calculate the numbers here and write out an initial \(\theta\) file instead.
The approx for \(\tau_f\) is known to be pretty rough, so after an initial run showed that the final \(\tau\) was about 2.6 MPa too high, I adjusted the \(\tau_f\) in the formula for \(\mu_0\) above by this amount (in code block below).
For TPV14 (and TPV15) we have that \(\tau_p = \mu_s \sigma = 0.677 \cdot 120 \mbox{ MPa} = 81.24 \mbox{ MPa}\) and \(\tau_f = \mu_d \sigma = 0.525 \cdot 120 \mbox{ MPa} = 63 \mbox{ MPa}\).
fricTheta0 = function(flt, rundir, mus=0.677, mud=0.525, sigma=120,
VEQ=1.22, Vstar=1e-6, Dc=1e-5, a=0.01, fa=0.1, b=0.015,
dtauf = -2.6,
ynuc=8000, znuc=-7500, wnuc=1500, segnuc=1) {
taup = mus * sigma[1]
tauf = mud * sigma[1]
aprime = fa*a
mu0 = (tauf + dtauf)/sigma[1] + (b-a)*log(VEQ/Vstar)
theta0 = (Dc/Vstar)^((b-a)/(b-a+aprime)) * ((l*b*sigma[1])/(lameMu*VEQ))^((aprime)/(b-a+aprime)) *
exp( (taup/sigma[1] - mu0)/(b - a + aprime) )
theta.nuc = RStheta(Vstar, mu0, taup, sigma[1], a, b, VEQ, Dc)
theta = rep(theta0, flt$np)
nuc = which(flt$iseg == segnuc & abs(flt$z - znuc) <= wnuc & abs(flt$y - ynuc) <= wnuc)
theta[nuc] = theta.nuc
write(theta, file=file.path(rundir, "initTheta.txt"))
return(mu0)
}
buildRSQSimInputFile =
function(rundir, run, a, b, Dc, mu0, ddotStar, VEQ, lameLambda, lameMu) {
unlink("eqs..out")
system2("/Users/dinger/RSQSim.github.master/runRSQSim", stdout=NULL, stderr=NULL)
eqs = readEqs("eqs..out", saveRData=FALSE)
unlink(c("eqs..out", "eqs..out.RData", ".Screen.Time.out", ".times_saved.out"))
params = eqs$params
rm(eqs)
params$A_1 = a
params$B_1 = b
params$Dc_1 = Dc
params$mu0_1 = mu0
params$ddotStar_1 = Vstar
params$alpha_1 = 0.25
params$ddotEQ_1 = VEQ
params$lameLambda = lameLambda
params$lameMu = lameMu
params$nEq = 1
params$faultFname = "faults.flt"
params$outFnameInfix = paste0("TPV14.run", run)
params$writeTau = 2
params$writeSigma = 2
params$initTauFname = "initTau.txt"
params$initSigmaFname = "initSigma.txt"
params$initThetaFname = "initTheta.txt"
writeRSQSimInputFile(params, file.path(rundir, paste0("TPV14.run", run, ".in")))
}
Note that with 1 km elements and set up as below, rupture would not propagate away from nucleation region. Currently using 400 m elements.
run = 1
rundir = file.path("TPV14", paste0("Run", run))
system(paste("mkdir -p", rundir))
l = 400
setoff = 0
ynuc = 8000
znuc = -7500
wnuc = 1500
segnuc = 1
taunuc = 81.6
mus = 0.677
mud = 0.525
initSigma = c(120, 120)
initTau = c(70, 70)
VEQ = 1.22 # note: usual default is 1.0, shear impedance in this case says 1.57
Vstar = 1e-6
Dc = 1e-5
a = 0.01
fa = 0.1
b = 0.015
dtauf = -2.6
eqs = NULL
flt = buildFault(l=l, setoff=setoff, rundir=rundir)
initStresses(flt, initSigma, initTau,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc, taunuc=taunuc)
mu0 = fricTheta0(flt, rundir, mus=mus, mud=mud, sigma=initSigma,
VEQ=VEQ, Vstar=Vstar, Dc=Dc, a=a, fa=fa, b=b,
dtauf = dtauf,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc)
buildRSQSimInputFile(rundir=rundir, run=run, a=a, b=b, Dc=Dc, mu0=mu0,
ddotStar=ddotStar, VEQ=VEQ,
lameLambda=lameLambda, lameMu=lameMu)
## [1] "Warning: /Users/dinger/Proposals/2019SCEC-branch/SCEC_CVS/faults.in doesn't exist"
here = getwd()
setwd(rundir)
unlink("eqs*out")
system2("mpirun",
args=c("-np", 4, "/Users/dinger/RSQSim.github.master/runRSQSim",
paste0("TPV14.run", run, ".in")),
stderr=paste0("run", run, ".stderr"),
stdout=paste0("run", run, ".stdout"))
eqs[[run]] = readEqs(paste0("eqs.TPV14.run", run, ".out"))
setwd(here)
hist(eqs[[1]]$taupList, breaks="FD")
c(mean(eqs[[1]]$taupList), mus*initSigma[1])
## [1] 81.18439 81.24000
tauf.rsq = readFaultVal(eqs[[1]], val="tau", writeVal="final")$val[1,]
hist(tauf.rsq, breaks="FD")
c(mean(tauf.rsq), mud*initSigma[1])
## [1] 62.84159 63.00000
run = 2
rundir = file.path("TPV14", paste0("Run", run))
system(paste("mkdir -p", rundir))
setoff = 100
flt = buildFault(l=l, setoff=setoff, rundir=rundir)
initStresses(flt, initSigma, initTau,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc, taunuc=taunuc)
mu0 = fricTheta0(flt, rundir, mus=mus, mud=mud, sigma=initSigma,
VEQ=VEQ, Vstar=Vstar, Dc=Dc, a=a, fa=fa, b=b,
dtauf = dtauf,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc)
buildRSQSimInputFile(rundir=rundir, run=run, a=a, b=b, Dc=Dc, mu0=mu0,
ddotStar=ddotStar, VEQ=VEQ,
lameLambda=lameLambda, lameMu=lameMu)
## [1] "Warning: /Users/dinger/Proposals/2019SCEC-branch/SCEC_CVS/faults.in doesn't exist"
here = getwd()
setwd(rundir)
unlink("eqs*out")
system2("mpirun",
args=c("-np", 4, "/Users/dinger/RSQSim.github.master/runRSQSim",
paste0("TPV14.run", run, ".in")),
stderr=paste0("run", run, ".stderr"),
stdout=paste0("run", run, ".stdout"))
eqs[[run]] = readEqs(paste0("eqs.TPV14.run", run, ".out"))
setwd(here)
run = 3
rundir = file.path("TPV14", paste0("Run", run))
system(paste("mkdir -p", rundir))
setoff = 200
flt = buildFault(l=l, setoff=setoff, rundir=rundir)
initStresses(flt, initSigma, initTau,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc, taunuc=taunuc)
mu0 = fricTheta0(flt, rundir, mus=mus, mud=mud, sigma=initSigma,
VEQ=VEQ, Vstar=Vstar, Dc=Dc, a=a, fa=fa, b=b,
dtauf = dtauf,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc)
buildRSQSimInputFile(rundir=rundir, run=run, a=a, b=b, Dc=Dc, mu0=mu0,
ddotStar=ddotStar, VEQ=VEQ,
lameLambda=lameLambda, lameMu=lameMu)
## [1] "Warning: /Users/dinger/Proposals/2019SCEC-branch/SCEC_CVS/faults.in doesn't exist"
here = getwd()
setwd(rundir)
unlink("eqs*out")
system2("mpirun",
args=c("-np", 4, "/Users/dinger/RSQSim.github.master/runRSQSim",
paste0("TPV14.run", run, ".in")),
stderr=paste0("run", run, ".stderr"),
stdout=paste0("run", run, ".stdout"))
eqs[[run]] = readEqs(paste0("eqs.TPV14.run", run, ".out"))
setwd(here)
run = 4
rundir = file.path("TPV14", paste0("Run", run))
system(paste("mkdir -p", rundir))
setoff = 300
flt = buildFault(l=l, setoff=setoff, rundir=rundir)
initStresses(flt, initSigma, initTau,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc, taunuc=taunuc)
mu0 = fricTheta0(flt, rundir, mus=mus, mud=mud, sigma=initSigma,
VEQ=VEQ, Vstar=Vstar, Dc=Dc, a=a, fa=fa, b=b,
dtauf = dtauf,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc)
buildRSQSimInputFile(rundir=rundir, run=run, a=a, b=b, Dc=Dc, mu0=mu0,
ddotStar=ddotStar, VEQ=VEQ,
lameLambda=lameLambda, lameMu=lameMu)
## [1] "Warning: /Users/dinger/Proposals/2019SCEC-branch/SCEC_CVS/faults.in doesn't exist"
here = getwd()
setwd(rundir)
unlink("eqs*out")
system2("mpirun",
args=c("-np", 4, "/Users/dinger/RSQSim.github.master/runRSQSim",
paste0("TPV14.run", run, ".in")),
stderr=paste0("run", run, ".stderr"),
stdout=paste0("run", run, ".stdout"))
eqs[[run]] = readEqs(paste0("eqs.TPV14.run", run, ".out"))
setwd(here)
Only continued to this offset because this is one of our element sizes, which what TPV14 did - but with smaller elements
run = 5
rundir = file.path("TPV14", paste0("Run", run))
system(paste("mkdir -p", rundir))
setoff = 400
flt = buildFault(l=l, setoff=setoff, rundir=rundir)
initStresses(flt, initSigma, initTau,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc, taunuc=taunuc)
mu0 = fricTheta0(flt, rundir, mus=mus, mud=mud, sigma=initSigma,
VEQ=VEQ, Vstar=Vstar, Dc=Dc, a=a, fa=fa, b=b,
dtauf = dtauf,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc)
buildRSQSimInputFile(rundir=rundir, run=run, a=a, b=b, Dc=Dc, mu0=mu0,
ddotStar=ddotStar, VEQ=VEQ,
lameLambda=lameLambda, lameMu=lameMu)
## [1] "Warning: /Users/dinger/Proposals/2019SCEC-branch/SCEC_CVS/faults.in doesn't exist"
here = getwd()
setwd(rundir)
unlink("eqs*out")
system2("mpirun",
args=c("-np", 4, "/Users/dinger/RSQSim.github.master/runRSQSim",
paste0("TPV14.run", run, ".in")),
stderr=paste0("run", run, ".stderr"),
stdout=paste0("run", run, ".stdout"))
eqs[[run]] = readEqs(paste0("eqs.TPV14.run", run, ".out"))
setwd(here)
run = 6
rundir = file.path("TPV14", paste0("Run", run))
system(paste("mkdir -p", rundir))
l = 200
setoff = 0
VEQ = 1.22 # note: usual default is 1.0, shear impedance in this case says 1.57
b = 0.015
dtauf = -2.6
flt = buildFault(l=l, setoff=setoff, rundir=rundir)
initStresses(flt, initSigma, initTau,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc, taunuc=taunuc)
mu0 = fricTheta0(flt, rundir, mus=mus, mud=mud, sigma=initSigma,
VEQ=VEQ, Vstar=Vstar, Dc=Dc, a=a, fa=fa, b=b,
dtauf = dtauf,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc)
buildRSQSimInputFile(rundir=rundir, run=run, a=a, b=b, Dc=Dc, mu0=mu0,
ddotStar=ddotStar, VEQ=VEQ,
lameLambda=lameLambda, lameMu=lameMu)
## [1] "Warning: /Users/dinger/Proposals/2019SCEC-branch/SCEC_CVS/faults.in doesn't exist"
here = getwd()
setwd(rundir)
unlink("eqs*out")
system2("mpirun",
args=c("-np", 4, "/Users/dinger/RSQSim.github.master/runRSQSim",
paste0("TPV14.run", run, ".in")),
stderr=paste0("run", run, ".stderr"),
stdout=paste0("run", run, ".stdout"))
eqs[[run]] = readEqs(paste0("eqs.TPV14.run", run, ".out"))
setwd(here)
hist(eqs[[run]]$taupList, breaks="FD")
c(mean(eqs[[run]]$taupList), mus*initSigma[1])
## [1] 81.18197 81.24000
tauf.rsq = readFaultVal(eqs[[run]], val="tau", writeVal="final")$val[1,]
hist(tauf.rsq, breaks="FD")
c(mean(tauf.rsq), mud*initSigma[1])
## [1] 62.73737 63.00000
run = 7
rundir = file.path("TPV14", paste0("Run", run))
system(paste("mkdir -p", rundir))
setoff = 100
flt = buildFault(l=l, setoff=setoff, rundir=rundir)
initStresses(flt, initSigma, initTau,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc, taunuc=taunuc)
mu0 = fricTheta0(flt, rundir, mus=mus, mud=mud, sigma=initSigma,
VEQ=VEQ, Vstar=Vstar, Dc=Dc, a=a, fa=fa, b=b,
dtauf = dtauf,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc)
buildRSQSimInputFile(rundir=rundir, run=run, a=a, b=b, Dc=Dc, mu0=mu0,
ddotStar=ddotStar, VEQ=VEQ,
lameLambda=lameLambda, lameMu=lameMu)
## [1] "Warning: /Users/dinger/Proposals/2019SCEC-branch/SCEC_CVS/faults.in doesn't exist"
here = getwd()
setwd(rundir)
unlink("eqs*out")
system2("mpirun",
args=c("-np", 4, "/Users/dinger/RSQSim.github.master/runRSQSim",
paste0("TPV14.run", run, ".in")),
stderr=paste0("run", run, ".stderr"),
stdout=paste0("run", run, ".stdout"))
eqs[[run]] = readEqs(paste0("eqs.TPV14.run", run, ".out"))
setwd(here)
run = 8
rundir = file.path("TPV14", paste0("Run", run))
system(paste("mkdir -p", rundir))
setoff = 200
flt = buildFault(l=l, setoff=setoff, rundir=rundir)
initStresses(flt, initSigma, initTau,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc, taunuc=taunuc)
mu0 = fricTheta0(flt, rundir, mus=mus, mud=mud, sigma=initSigma,
VEQ=VEQ, Vstar=Vstar, Dc=Dc, a=a, fa=fa, b=b,
dtauf = dtauf,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc)
buildRSQSimInputFile(rundir=rundir, run=run, a=a, b=b, Dc=Dc, mu0=mu0,
ddotStar=ddotStar, VEQ=VEQ,
lameLambda=lameLambda, lameMu=lameMu)
## [1] "Warning: /Users/dinger/Proposals/2019SCEC-branch/SCEC_CVS/faults.in doesn't exist"
here = getwd()
setwd(rundir)
unlink("eqs*out")
system2("mpirun",
args=c("-np", 4, "/Users/dinger/RSQSim.github.master/runRSQSim",
paste0("TPV14.run", run, ".in")),
stderr=paste0("run", run, ".stderr"),
stdout=paste0("run", run, ".stdout"))
eqs[[run]] = readEqs(paste0("eqs.TPV14.run", run, ".out"))
setwd(here)
run = 9
rundir = file.path("TPV14", paste0("Run", run))
system(paste("mkdir -p", rundir))
setoff = 300
flt = buildFault(l=l, setoff=setoff, rundir=rundir)
initStresses(flt, initSigma, initTau,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc, taunuc=taunuc)
mu0 = fricTheta0(flt, rundir, mus=mus, mud=mud, sigma=initSigma,
VEQ=VEQ, Vstar=Vstar, Dc=Dc, a=a, fa=fa, b=b,
dtauf = dtauf,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc)
buildRSQSimInputFile(rundir=rundir, run=run, a=a, b=b, Dc=Dc, mu0=mu0,
ddotStar=ddotStar, VEQ=VEQ,
lameLambda=lameLambda, lameMu=lameMu)
## [1] "Warning: /Users/dinger/Proposals/2019SCEC-branch/SCEC_CVS/faults.in doesn't exist"
here = getwd()
setwd(rundir)
unlink("eqs*out")
system2("mpirun",
args=c("-np", 4, "/Users/dinger/RSQSim.github.master/runRSQSim",
paste0("TPV14.run", run, ".in")),
stderr=paste0("run", run, ".stderr"),
stdout=paste0("run", run, ".stdout"))
eqs[[run]] = readEqs(paste0("eqs.TPV14.run", run, ".out"))
setwd(here)
run = 10
rundir = file.path("TPV14", paste0("Run", run))
system(paste("mkdir -p", rundir))
setoff = 400
flt = buildFault(l=l, setoff=setoff, rundir=rundir)
initStresses(flt, initSigma, initTau,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc, taunuc=taunuc)
mu0 = fricTheta0(flt, rundir, mus=mus, mud=mud, sigma=initSigma,
VEQ=VEQ, Vstar=Vstar, Dc=Dc, a=a, fa=fa, b=b,
dtauf = dtauf,
ynuc=ynuc, znuc=znuc, wnuc=wnuc, segnuc=segnuc)
buildRSQSimInputFile(rundir=rundir, run=run, a=a, b=b, Dc=Dc, mu0=mu0,
ddotStar=ddotStar, VEQ=VEQ,
lameLambda=lameLambda, lameMu=lameMu)
## [1] "Warning: /Users/dinger/Proposals/2019SCEC-branch/SCEC_CVS/faults.in doesn't exist"
here = getwd()
setwd(rundir)
unlink("eqs*out")
system2("mpirun",
args=c("-np", 4, "/Users/dinger/RSQSim.github.master/runRSQSim",
paste0("TPV14.run", run, ".in")),
stderr=paste0("run", run, ".stderr"),
stdout=paste0("run", run, ".stdout"))
eqs[[run]] = readEqs(paste0("eqs.TPV14.run", run, ".out"))
setwd(here)